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ABSTRACT 

We analyze the gas accretion flow through a planet-produced gap in a protoplanetary disk. We adopt the 
alpha disk model and ignore effects of planetary migration. We develop a semi-analytic, one-dimensional 
model that accounts for the effects of the planet as a mass sink and also carry out two-dimensional 
hydrodynamical simulations of a planet embedded in a disk. The predictions of the mass flow rate through 
the gap based on the semi-analytic model generally agree with the hydrodynamical simulations at the 25% 
level. Through these models, we are able to explore steady state disk structures and over large spatial 
ranges. The presence of an accreting ~ 1 Mj planet significantly lowers the density of the disk within a 
region of several times the planet’s orbital radius. The mass flow rate across the gap (and onto the central 
star) is typically 10% to 25% of the mass accretion rate outside the orbit of the planet, for planet-to-star 
mass ratios that range from 5 x 10~ 5 to 1 x 10 -3 . 

Subject headings: accretion, accretion discs - hydrodynamics — planets and satellites: general 


1. Introduction 

The presence of a ~ 1 Mj planet influences the structure of a 
circumstellar disk. The most obvious effect is the opening of a 
tidally produced gap (Lin & Papaloizou 1993). Recent stud¬ 
ies of young stars such as CoKu Tau/4 (D’Alessio et al. 2005), 
TW Hya (Calvet et al. 2002), KH 15D (Herbst et al. 2002), 
HR 4796A (Schneider 1999), and HD 141569A (Clampin et 
al. 2003) suggest the presence of inner holes in circumstellar 
disks. One possible mechanism for producing such a hole is 
the tidal barrier created by a planet. A planet’s tidal forces 
can, in principle, prevent material from accreting across the 
orbit of the planet. The unreplenished material interior to the 
planet’s orbit accretes onto the central star, thereby produc¬ 
ing a hole. However, alternative explanations for apparent 
holes have been offered for some cases, such as dust segrega¬ 
tion (Takeuchi & Artymowicz 2001). 

This paper deals with the determination of the extent to 
which the disk material exterior to a planet’s orbit can re¬ 
plenish the interior material. A planet of ~ 1 Mj opens a 
gap in a disk having typical parameters. But the gap is not 
clean and is accompanied by accretion onto the planet and 


1 To appear in The Astrophysical Journal sometime in 2006. 


accretion flow across the gap (Artymowicz & Lubow 1996; 
Bryden et al. 1999; Kley 1999; Lubow, Seibert, & Artymow¬ 
icz 1999). Several other past studies of planet-disk interac¬ 
tions have concentrated on analyzing torques and accretion 
flows onto planets. They also have sometimes reported find¬ 
ing mass flow through gaps (e.g., Winters, Balbus, & Hawley 
2003), while other studies have reported little inflow (e.g., 
Bate et al. 2003). 

The extent of the flow through the gap depends on the 
planet mass and disk properties. In the study of a 1 Mj planet 
that orbits a 1 Mg star by Lubow, Seibert, & Artymowicz 
(1999), very little mass was found to accrete interior to the 
planet’s orbit, not surprisingly, with nonaccreting boundary 
conditions at the inner boundary located at 0.3 times the 
orbital radius of the planet. On the other hand, if the interior 
disk was initially strongly depleted of material, then there was 
a substantial accretion flow that built up the interior disk. In 
this configuration, the accretion rate into the interior disk is 
comparable to the rate onto the planet. A similar density 
profile was considered recently by Quillen et al. (2004) to 
model CoKu Tau/4. Consequently, the accretion rate onto 
the central star depends on the density distribution and the 
inner boundary conditions. 

Even low mass planets which do not open gaps can affect 
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the global structure of a circumstellar disk. The previous 
studies of flow across a gap are largely based on multidimen¬ 
sional hydrodynamical simulations. They are limited by their 
dynamic range in space, which typically involve a region that 
ranges from 0.4 to 6 times the orbit radius of the planet. 
A smaller inner boundary increases the computational over¬ 
head, due the need to take shorter timesteps because of the 
higher rotational velocity. Simulations are limited in time to 
at most a few thousand orbits of the planet, which is less than 
the viscous evolution timescale of a typical protostellar disk. 
To overcome these limitations and provide some physical in¬ 
sight, we develop in Section 2 a semi-analytic model for the 
steady state flow across the gap. This model depends on a 
key parameter, the accretion efficiency, which we determine 
by two-dimensional simulations described in Section 3. The 
semi-analytic model predicts the rate of accretion past the 
planet and the steady state density profile. The predictions 
are compared with results of numerical simulations. Section 4 
contains the conclusions. 

2. Semi-analytic Model 

We develop a one-dimensional, semi-analytic model for a cir¬ 
cumstellar disk containing a planet. We assume that the 
planet lies on a fixed circular orbit and neglect effects of mi¬ 
gration. In this model, we separate the gap region from the 
main disk region. We assume that the tidal torques on the 
disk exerted by the planet are confined to the gap region. 
In reality, waves generated by the planet propagate into the 
main disk and exert additional torques, which we ignore. But, 
shocks and other processes limit the extent of wave propaga¬ 
tion (Bate et al. 2002). 

We adopt the alpha disk model for the disk turbulence 
throughout. If the disk turbulence is due to magnetic fields, 
then the presence of a planet can affect the nature of the 
turbulence (Nelson & Papaloizou 2003; Winters, Balbus, & 
Hawley 2003). We do not account for such effects. 

We model the effect of the planet on the main disk as 
a mass sink that lies on the planet’s orbit. The mass sink 
description utilizes the accretion efficiency parameter 

£ =^k- (1) 

where M p is the accretion rate onto the planet and v p is the 
turbulent kinematic viscosity. Density E p is the disk density 
at the location of the planet, based on a smooth continuation 
(interpolation) of the density profile outside the gap to the 
location of the planet (Lubow, Seibert, & Artymowicz 1999). 
For a narrow gap, this density is the density just outside the 
gap. The denominator is recognized as the standard from for 
the steady state rate of accretion through a viscous disk, far 
from the disk inner boundary (Lynden-Bell & Pringle 1974). 
Since the accretion rates scale with v , we expect and find 
that £ is fairly independent of v (and E p ) for fixed H/r, as 
is consistent with the results of Kley (1999). The efficiency 
depends mainly on the planet mass M p . One worrisome result 
of the numerical simulations is that £ is found to be greater 


than unity. One might think that this suggests the planet 
is accreting more mass than the disk can supply in a steady 
state. We show later that steady state accretion flows can 
and usually do have £ greater than unity. 

The density structure of a circumstellar accretion disk 
evolves in time from some initial state. The initial disk state 
of a protostellar disk is not known. Disks tend to evolve 
towards a steady state in which the accretion rate is indepen¬ 
dent of radius (Lynden-Bell & Pringle 1974; Pringle 1981). 
The disk evolutionary timescale in the vicinity of planets 
that lie within a 10 AU of the central star is shorter than 
the global evolutionary timescale for a 100 AU disk, and so a 
steady state will be even more likely achieved on these smaller 
scales. Consequently, we seek steady state solutions for the 
disk structure. 


2.1. One-Dimensional Equations 


With the above model, the one-dimensional equations of mass 
conservation and azimuthal force balance for a disk with ra¬ 
dial velocity u, angular velocity D, and surface density E in 
a steady state (dE/dt = 0 ,du/dt = 0,dv/dt = 0) Keplerian 
disk are 


1 d (rEu) 
r dr 


3uE£ 
2 r 


8(r 


r p) 


( 2 ) 


and 

r 2 £lEu = _ d ^l lQr 1 + 2rSA(r), (3) 

dr 

where p = vE with kinematic turbulent viscosity v, and 8 is 
the Dirac delta function. A(r) is the torque density per unit 
mass produced by the tidal field of the planet 


A(r) = Sign(r 


fq 2 GM* 
2 r 



(4) 


where / is a constant of order unity and A p is the maximum 
of H and |r — r p \. These equations are identical to those in 
Lin & Papaloizou (1986), except that we have added to sink 
term in the equation of continuity (2). 

The determination of the accretion efficiency £ in equa¬ 
tion (2) requires a model for the capture of gas by the planet. 
We shall rely on two-dimensional simulations to determine 
its value. We use the one-dimensional model to provide ini¬ 
tial conditions for the two-dimensional simulations. To com¬ 
pose the one-dimensional model, we divide space into two re¬ 
gions: the gap region and the main disk (nongap) region. We 
combine the solutions in the two regions of space to obtain 
approximate steady state global solutions. In the process, 
we self-consistently determine the accretion efficiency, which 
links the solutions in these two regions. 


2.2. Gap Region 

From previous two and three-dimensional studies discussed in 
the Introduction, we know that the gap region is complicated 
by the presence of nonKeplerian flow and shocks. We do 
not expect that a one-dimensional model would work well in 
describing the detailed gas flow there. But, we estimate the 
density structure E g in the gap region by assuming d/dr 
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l/r (WKB approximation) and 1/A P 1/r in equation (3) 
to obtain 


E 9 (r) = E p exp 


fq 2 r 2 p n p 

b 3 

CO 

9 v p 

>1 
^ 1 

1 _ 


For a narrow gap, it follows that 

fu(xp) = [l e {Xp), 


( 11 ) 


(5) 


where E p is the density that would occur at r = r v in the ab¬ 
sence of a gap, as defined in equation (1). We generally take 
constant / equal to 2. We have used initial density profiles like 
this in our previous work because they provide good initial 
conditions for planet-disk simulations. They undergo moder¬ 
ate changes and appear to reach a dynamically steady state 
in the gap after a few hundred planetary orbits. The achieve¬ 
ment of a dynamically steady state is expedited through the 
development of shocks. 

On longer, viscous timescales, the gap region could un¬ 
dergo further structural changes. They arise in part because 
of the interaction of the gap region with the main disk region. 
The coupling between these two regions occurs through the 
accretion efficiency parameter £, which is determined in the 
gap region by the outcome of the simulations. For a consistent 
solution, we require that the value of the accretion efficiency 
parameter match in the two regions, as described in the next 
subsection. 

2.3. Main Disk Region 

Outside the gap region, we ignore the presence of waves which 
we assume damp most of their energy near the planet and 
neglect the planetary torque A. The planetary accretion af¬ 
fects the large scale accretion flow as a mass sink. We expect 
that main disk region is described well by the one-dimensional 
model. 

From equation (2), we have that 


since fii(x p ) and /i e (x p ) are nearly equal to the fj,- values on 
both sides of the planet just outside the gap. We combine 
equations (8), (9), and (10) to obtain 


Aq )x) 


x p (x — x t ) 
x(x p — x t ) 


and 


/te (a;) = Atp 


I + ('-?)( f + 


X* 


( 12 ) 


(13) 


where fi p is the value of a t smoothly extended from the main 
disk to the planet, excluding the gap, fi p = Hi(x p ) = fj, e (x p ). 

Notice that for £ = 0 (no planet) and x t = 0, the solution 
implies that fi(x) is constant, which is equivalent to the con¬ 
dition that M = 3tyvT, is independent of r in a steady state, 
as discussed above. 

Far outside the orbit of the planet at some x 0 x p , the 
disk is hardly aware of the planet. We have from equations (3) 
and (13) that 


Me = 37r/Xo = 37r/Xp £ + 


(14) 


Exterior accretion rate M e is independent of x for all x > 


x p . The accretion rate onto the planet is given by 
M p = 37ta i p £, 


(15) 


and the accretion rate interior to the planet’s orbit (and onto 
the star) is then 


(16) 


Mi = —2nrUiY,i 

(6) 

II 

1 

■<2 

II 

37r fjjp Xp 

Me = —27rru e E e 

(7) 

Xp — X* 

M e = Mi + 3ni'pY,p£, 

(8) 

which is constant for x < x p . 



where M; and M e denote the mass accretion rate respectively 
interior and exterior to the orbit of the planet, and u p and E p 
are as used in equation (1). 

We solve equation (3) subject to the zero-stress boundary 
condition at the inner edge to account for the effects of a 
boundary layer, where the disk meets the central star, as 
described by Lynden-Boll & Pringle (1974). To model this 
effect, one typically applies a zero density boundary condition 
at the inner edge (e.g., Pringle, Verbunt & Wade 1986), so 
that A t ( r *) = 0, where r, is the location of the stellar radius. 
The solution is given by 


The mass flow rate ratio of accretion past the planet to 
accretion onto the planet is then 


Mi. 

M 


P (1 i / ft/r p )£ 


(17) 


The ratio of the accretion rate interior to the orbit of the 
planet to the rate exterior is given by 


Mi_ 

Me 


1 


i + (i - V r */ r p)£ 


(18) 


3niu(x) = Mil 1 — 


and 


3n/j. e (x) = Mo H- 


-) 

(9) 

X* J 


C 


X 

(10) 


where x = t/r, fj. p = fj,(x p ), and /Xi (A t e) is the value of [i inte¬ 
rior (exterior) to the orbit of the planet, and C is a constant 
of integration. 


Ratios (17) and (18) are independent of fi p . 

The star and planet compete for accretion flow. The above 
accretion rate ratios indicate that the closer the stellar sur¬ 
face comes to the planet, the more the star accretes and the 
less the planet accretes. This prediction is confirmed in the 
numerical simulations discussed later (model b versus model 
g). This result is a consequence of the influence of the star 
in diverting flow from the planet. For a given kinematic vis¬ 
cosity v(r ) and accretion efficiency £, we can determine the 
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disk density profile outside the gap £(r) = p(r)/v(r ) from 
equations (12) and (13). The accretion efficiency depends on 
the flow details within the gap, which we determine in Sec¬ 
tion 3 based on two-dimensional hydrodynamical simulations. 
Equation (17) determines the ratio of the accretion rate onto 
the central star to the accretion rate onto the planet. This ra¬ 
tio can be compared against results of numerical simulations. 

For high mass planets (typically several Jupiter masses), 
the analytic model breaks down because gravitational torques 
cannot be ignored in the main disk (outside the gap). This 
breakdown can be seen in equation (18). For a high planet 
mass, the gap is very clean and the accretion efficiency £ is 
very small. In the limit of a (nonmigrating) high planet mass, 
the disk will evolve towards a decretion disk, rather than an 
accretion disk (Pringle 1981). A decretion disk has a very 
different description than the accretion disk description used 
here. In the high planet mass limit, equation (18) predicts 
that mass accretes uninhibitedly past the planet, which is in¬ 
correct. Instead, equation (18) applies to systems containing 
lower mass planets. In that case, £ is small when the planet 
mass is small and equation (18) provides a proper description. 

Notice that there is no problem with having £ greater than 
unity in a steady state. The accretion rate of material outside 
the orbit of the planet, M e , is greater than 3np p . The dy¬ 
namic viscosity and typically the surface density near the 
planet (and outside the gap) are reduced due to the presence 
of the planet, in accordance with equation (14). 

3. Numerical Hydrodynamical Simulations 
3.1. Description of Code 

We carried out a series of numerical disk simulations in two 
dimensions using the code described in D’Angelo, Henning, 
& Kley (2002); D’Angelo, Kley, & Henning (2003). The code 
calculates the time-evolution Navier-Stokes equations for a 
gaseous disk that or bits a central star. The disk is subject to 
gravitational forces from a planet that lies on a fixed circular 
orbit. The code allows for high resolution near the planet 
by using set of nested grids. We adopt a locally isother¬ 
mal equation of state with sound speed equal to the disk 
aspect ratio, H/r, times the local Keplerian velocity. We 
used a constant disk aspect ratio H/r = 0.05 throughout, 
implying that the temperature distribution scales as the in¬ 
verse of the distance from the disk axis. Unless otherwise 
stated, the kinematic viscosity at the radius of the planet 
is Up = 10~ 5 r p f2 p , which is equivalent to Shakura & Sun- 
yaev parameter a = 4 x 10~ 3 at the location of the planet. 
We consider a range of planet masses, having a mass ratio 
with the central star of q = 5 x 10 -5 to 2 x 10~ 3 . Be¬ 
low about q = 4 x 10 -4 with a lMg star, the size of the 
planet’s Hill sphere is smaller than the disk thickness and 
three-dimensional effects can be important (Bate et al. 2003). 
Due to computational resource limitations, the simulations 
were carried out in two dimensions. In low end of the planet 
mass range we consider, numerical estimates of the accre¬ 
tion rate in two and three dimensions differ by only 30% 
(D’Angelo, Kley, & Henning 2003). So we do not expect that 


three-dimensional effects will substantially alter our conclu¬ 
sions. 

The grid size used for these calculations consisted of a 
three-level nested-grid system with the basic level having 
374 x 422 grid zones in the radial and azimuthal direction, 
respectively. The first and the second subgrid level had 
113 x 103 and 133 x 123 grid zones, respectively. The highest 
resolution achieved in the gap region around the planet was 
Ar = 3.7 x 10 _3 r p . The resolution in the azimuthal direction 
was equal to A r/r p . 

The origin of the coordinate system is located on the star 
and the reference frame rotates about the disk axis at a rate 
equal to the angular velocity of the planet. The acceleration of 
the coordinate system origin relative to center-of-mass of the 
star-planet system is accounted for in the gravitational po¬ 
tential of the disk. As discussed in Section 3.3, we performed 
a convergence test on one of our models which indicates that 
this resolution is adequate. Self-gravity of the disk material 
is ignored. 

The smoothing length of the planet potential was chosen 
to be 0.1 times the planet’s Hill radius. Two and three- 
dimensional models with that smoothing length yield similar 
accretion rates around a Jupiter mass planet (D’Angelo, Kley, 
& Henning 2003). Some recent models show that variations 
of this smoothing by a factor of about 1.4 affects the accretion 
rate by about 1%. 

Accretion boundary conditions were generally applied to 
the planet. Such boundary conditions are currently thought 
to be appropriate to the range of planet masses considered in 
this paper because such planets are believed to undergo run¬ 
away gas accretion (e.g., Pollack et al. 1996). The accretion 
was simulated by removing mass within a radius of 0.2 Hill 
radii from the planet. We fix a removal time-scale, which is 
on the order of 0.1 orbital periods. Tests on how the accretion 
rate depends on these two parameters are given in D’Angelo, 
Henning, & Kley (2002). The tests show that the accretion 
rate does not depend sensitively on the details of the mass re¬ 
moval parameter values, provided the removal prevents mass 
build-up near the planet. The formulation of planetary mass 
accretion in Lubow, Seibert, & Artymowicz (1999) and Bate 
et al. (2003) is completely different from that adopted here. 
Nonetheless, they obtain accretion rates that are very simi¬ 
lar to those obtained by (D’Angelo, Henning, & Kley 2002; 
D’Angelo, Kley, & Henning 2003). Consequently, the evi¬ 
dence all suggests that the accretion rate is the rate that a 
planet would accrete under runaway accretion conditions. It 
is not an artifact of the mass removal process. 

The inner boundary was usually set at 0.4r p , while the 
outer boundary was located at 6 r v . At the inner boundary, 
outflow boundary conditions were applied. This means that 
if u r < 0 at the first active zone, this value is transferred 
to the zones off the active grid (ghost zones). If u r > 0 at 
the first active zone, then u r of the ghost cells is set to zero. 
At the outer boundary, inflow and outflow of material are 
permitted. The outflow is implemented as in Godon (1996, 
1997). During the course of the simulation, the values of all 
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physical variables are free to change at the inner and outer 
boundaries. They are not fixed by the initial analytic model. 


3.2. Density Determination 

For a fixed choice of v(r) and an initial guess of accretion ef¬ 
ficiency £, we determine the initial density profile from equa¬ 
tions (12) and (13) with a superimposed gap profile (5). This 
profile is given by 


E(r) 


M r ) s g( r ) 

v(r)T, p 


(19) 


where /x(r) refers to fM for r < r v and /r e for r > r p . 

We generally set r* according to the location of the inner 
boundary. The initial gap profile is determined by an ap¬ 
proximate balance of viscous torques with tidal torques, as 
described in Section 2.2. The value of £ in the simulations 
is determined by the accretion rate onto the planet through 
equation (15). We run the simulation for typically 700-800 
orbits at which point a nearly time-independent value of £ 
is determined (see Figure 1). However, this value in general 
disagrees with the value assumed in constructing the den¬ 
sity profile. We then guess another value of £, construct the 
corresponding density profile, and rerun the simulation. We 
continue iterating until the assumed value of £ used for the 
density profile matches the measured value from equation (15) 
to within 10%. This process typically requires 3 or 4 itera¬ 
tions. 

Although we find a nearly stationary value of £ (Figure 1), 
the determination of a strict steady state profile in the gap 
from arbitrary initial conditions would require running the 
code beyond a local viscous timescale of order 10 4 orbits, 
much longer than 700 orbits. However, the gap structure 
is also determined by nonlinear effects (shocks) that act on 
shorter timescales and are likely responsible for the near- 
stationarity of £. Furthermore, from Figure 1 we see that 
£ converges in time to the same value, independent of initial 
conditions. This gives us confidence that the value of £ we de¬ 
termine from simulations is the steady state value. Nonethe¬ 
less, we cannot prove that the value of £ that we determine 
is rigorously the steady state value without integrating over 
longer timescales. 


3.3. Description of Results 

The results of the simulations are contained in Table 1. 
The default parameters consist of a kinematic viscosity v = 
10“ 5 r 2 fl p which is taken to be independent of radius, to¬ 
gether with the other parameters described in Section 3.1. 
The first column in Table 1 contains the model label. The 
second column contains the planet mass ratio. The third 
column contains the variant, if any, on the default parame¬ 
ters. The fourth and fifth columns contain respectively the 
accretion efficiency £ applied to the final iteration and the 
accretion efficiency £ obtained from the final simulation. As 
discussed in Section 3.2, for each model we iterated on £ un¬ 
til the applied and simulated £ agreed to within 10%. The 
sixth column contains the ratio of the accretion rate past the 


planet (and onto the central star) to the accretion rate onto 
the planet as obtained from the simulation. The seventh col¬ 
umn contains the same quantity as the sixth, but obtained 
from equation (17) of the semi-analytic model. The eighth 
column contains a correction, described in Section 3.6, that 
is applied to the values in the seventh column, in order to 
account for the location of the inner boundary. 

Model b is the default q = lx 1CP 3 model. The next seven 
models are variants of model b having a different initial den¬ 
sity profiles (models c and d), twice the kinematic viscosity 
(model e), a viscosity that increases as r 1 ^ 2 with the same 
H/r = 0.05 (model f), a smaller inner boundary (model g), a 
lower temperature H/r = 0.03 disk (model h), and a planet 
that does not accrete gas (model i). 

A convergence test was performed on model b by using 
a degraded resolution throughout the computational domain. 
The resolution was reduced by 25%, in each direction, relative 
to that presented in Section 3.1. Both in terms of accretion 
efficiency and accretion rate past the planet, the outcomes 
agree within 2% with those in the Table. 



Orbits 

Fig. 1.— Plot of the averaged accretion efficiency (£), aver¬ 
aged over the previous 10 planet orbits, versus time in planet 
orbit periods for a disk containing a q = 1 x 10~ 3 planet. 
The models in the three cases have different initial density 
distributions in the gap region. The curves from lowest to 
highest correspond to models b, c, and d of Table 1, which 
have values of / = 2,1, and 0.5, respectively for the initial 
gap density distribution of equation (5). 
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Table 1 

: Results from simulations (see 

Section 3.3) 




Model 



q 

variant 

Applied 

£ 

Simulated 

£ 

Simulated 

Mi/Mp 

Predicted 

Mi/Mp 

Corrected 

Mi/Mp 

a 

2 

X 

nr 3 


3 

2.9 

0.52 

0.91 


b 

1 

X 

icr 3 


8 

8.3 

0.38 

0.34 

0.13 

c 

1 

X 

icr 3 

/ = 1 

8 

8.3 

0.38 

0.34 

0.13 

d 

1 

X 

icr 3 

/ = 0.5 

8 

8.4 

0.39 

0.34 

0.13 

e 

1 

X 

10 -3 

2xi/ 

8 

8.8 

0.35 

0.34 

0.13 

f 

1 

X 

10 -3 

v oc r 1 ' 2 

8 

8.2 

0.40 

0.34 

0.13 

g 

1 

X 

10 -3 

ri = 0.22 r p 

9 

9.4 

0.14 

0.11 

0.11 

h 

1 

X 

icr 3 

H/r = 0.03 

6 

6.3 

0.29 

0.46 


i 

1 

X 

icr 3 

II 

o 

0 

0 




j 

5 

X 

icr 4 


12 

12.4 

0.35 

0.23 

0.08 

k 

1 

X 

icr 4 


7 

7.3 

0.50 

0.39 

0.14 

1 

5 

X 

10 -5 


3 

3.2 

1.10 

0.91 

0.33 


3.4. Accretion Efficiency 

As can be seen from the Table, the accretion efficiency is 
largely a function of planet mass and is insensitive to the ini¬ 
tial density profile (models c and d). It is also insensitive to 
the viscosity and its variation in r, and inner boundary loca¬ 
tion which are varied in models e through g, relative to model 
b. The efficiency varies from a low value of 3 for the highest 
and lowest mass planets considered, to a high value of 12 at 
an intermediate mass. This behavior is easily understood. At 
the high mass end, q = 2 x 10~ 3 (model a), the tidal field of 
the planet is so strong that material is inhibited from accret¬ 
ing onto the planet, resulting in a lower £ value relative to 
model b of q = 1 x 10~ 3 . 

However, the accretion across the gap is somewhat less 
effected, since Mi/M p is larger for q = 2 x 10“ 3 than for 
1 x 10~ 3 . While at the low mass end (model 1), the planet’s 
gravitational field has less ability to attract and accrete mat¬ 
ter. For comparison, the Jupiter-mass models (q = 1 x 10~ 3 ) 
in D’Angelo, Kley, & Henning (2003) and Bate et al. (2003) 
have an efficiency of about 4. As discussed in the Introduc¬ 
tion, a planet’s accretion rate is influenced by the density 
structure outside the gap (r > r p ), which in the present 
steady state case is different from that used in the above 
papers. 

3.5. Accretion Past Planet 

The last three columns in Table 1 describe the rate at which 
material flows past the orbit of the planet as a fraction of 
the accretion rate onto the planet. The values predicted by 
equation (17) generally agree with the values obtained by the 
simulations at the 25% level. 

We believe model a (q = 2 x 10 -3 ) is in the high mass 
regime where the semi-analytic model breaks down because 


gravitational torques in the main disk cannot be ignored, as 
discussed in Section 2.3 for decretion disks. The disk in this 
case is not a true decretion disk because there is some ac¬ 
cretion through the gap. However, the accretion efficiency is 
markedly lower than that in model b which has q = 1 x 10 -3 . 
We regard model a as a transitional disk between accretion 
and decretion. We believe this explains the discrepancy be¬ 
tween the simulated and predicted accretion ratios. 

For model b (q = 1 x 1CP 3 ) the simulated and predicted 
accretion rate ratios are in very good agreement, suggesting 
that the disk with such a planet behaves as an accretion disk. 
Agreement is also very good for model e in which the disk vis¬ 
cosity has been doubled and model f, in which the v increases 
with radius. 

Model g has a smaller inner boundary located at 0.22r p , 
instead of the default 0.4r p , and used an initial density profile 
that has r* = 0. It attempts to provide more realistic coverage 
of the interior disk, but its computational domain does not 
cover the assumed profile extent, as do all the other models 
(which have r* = 0.4r p ). Both the simulated and predicted 
flow rates indicate that the flow rate past the planet decreases 
significantly for a smaller inner boundary radius. 

Model h has a colder disk and consequently the tidal effects 
of the planet on the disk which open the gap are relatively 
stronger than pressure and viscosity. This case is intermediate 
between models a and b, with accretion efficiency lower than 
for model b. The simulated and predicted mass flow ratios 
in columns 6 and 7 do not agree well, probably for similar 
reasons given for model a. The flow onto the planet and 
across the gap are both reduced by about 25% relative to 
model b. 

Model i has nonaccreting boundary conditions onto a 
planet having q = lx ICC 3 . Consequently, we have applied an 
initial model having £ = 0 to the simulation. The accretion 
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rate obtained from the simulation agrees with the predicated 
steady state accretion rate, equation (16), to within 5%. The 
accretion rate onto the star should be unaffected by the pres¬ 
ence of the planet. 

The agreement on flow rates for model j is the least sat¬ 
isfactory of all models having H/r = 0.05 and q < lx 10~ 3 . 
It is not clear why this is the case. It may have to do with 
the high accretion efficiency, which is also observed in three- 
dimensional calculations (D’Angelo, Kley, & Henning 2003). 
Fairly good agreement (~ 20%) between the simulated and 
predicted accretion rates is found for the lowest mass models 
k and 1. 

3.6. Correction Due To Inner Boundary 

For planets whose orbital radii are of order AU, the inner 
boundary radius is relatively small, r* < 0.01r p . It is difficult 
to directly simulate a disk with a small inner boundary, due 
to the shortness of timesteps required. We estimate a cor¬ 
rection factor for the boundary location to be applied to the 
accretion rates. We saw in the previous subsection that the 
accretion efficiency £ depends mainly on planet mass and is 
insensitive to the inner boundary location. Using the default 
inner boundary location of 0.4r p and applying equation (17), 
we suggest that a correction factor to the flow rate ratio for a 
small radius inner boundary is about 1 — a/0.4 rs 0.37. This 
factor is applied to the predicted accretion rate ratios (except 
for models a and li for which the analytic model does not ap¬ 
ply, model i which does not allow accretion onto the planet, 
and model g which already has r* = 0). The result is in the 
eighth column of Table 1. 

The ratio of accretion rate onto the star to the rate onto 
the planet varies with mass ratio. The ratio is about 13% for 
a q = 1 x 10 -3 planet and 33% for a q = 5 x 1CD 5 planet. 
From equation (18) with r* = 0, the ratio of the accretion 
rate onto the star to the accretion rate outside the orbit of 
the planet is 11% for a q = 1 x 1CD 3 planet and 25% for a 
q = 5 x 1CT 5 planet. 

3.7. Density Profiles 

We examine the effect of a planet on the density distribution 
of a steady state viscous disk. Figure 2 shows the assumed 
density profile and evolved density profile for several models 
which are converged in £. For all the cases, the code uses 
the same boundary conditions, as discussed in Section 3.1. 
There are some differences in the profiles as a consequence 
of evolution, including some narrowing of the gap. Overall, 
however, the agreement between the analytic (initial) profile 
and the evolved one is very good. 

Global density profiles with a more realistic inner bound¬ 
ary location (r* = 0) were constructed for a q = 1 x 10~ 3 
planet using £ = 8 (see Table 1). A global analytic density 
distribution is obtained from equation (19) for some assumed 
viscosity function v(r). In addition, density profiles for cor¬ 
responding disks without a planet were also determined. The 
results are shown in Figure 3. These density distributions do 
not contain the evolutionary effects that in principle could 


be determined by numerical simulations. But as can be seen 
from Figure 2, these effects are modest, at least over 700-800 
orbits. Furthermore, it would be impossible to carry out nu¬ 
merical simulations with a realistically small inner boundary 
due to the shortness of timesteps required. The cases without 
planets correspond to steady disks that have the same viscos¬ 
ity distribution at all r and the same density distribution at 
r r p as the corresponding cases with planets. (For the 
left plot in Figure 3, the density match occurs at distances 
larger than shown). If an accreting planet were placed in a 
disk having the non-planet density distribution, indicated by 
dashed curves in Figure 3, then over time it would evolve to¬ 
wards the density distributions, indicated by the solid curves 
in Figure 3. 

The accreting planet effects the disk density distribution 
over several r p . For r r p , the fractional density difference 
between the nonplanet and planet cases varies slowly as ~ 
sjr p /r. Consequently, the planet’s effects on the disk extend 
over several times r p . Also plotted are the cases with an 
inner boundary of r* = 0.4r p , as was generally adopted in 
the numerical simulations. Notice that such cases agree well 
with the r« = 0 cases for r > r p , but poorly represent the 
density at small r. 

4. Discussion and Summary 

We have investigated two-dimensional steady state configura¬ 
tions of a protostellar disk containing a planet that undergoes 
mass accretion, using both an analytic model and numerical 
simulations. The planet is assumed to be of fixed mass and 
on a fixed circular orbit. Planet migration and mass gain can 
be neglected at late stages of evolution where the disk mass 
is less than the planet mass. The results also suggest the 
outcome for the more general case involving mass gain and 
migration. 

We speculate on the nonsteady effects of planet mass gain 
and migration. The mass doubling timescale for Jupiter un¬ 
dergoing runaway gas accretion from the minimum mass so¬ 
lar nebula is of order the local viscous timescale in our model 
(~ 10 s years). Consequently, the nonsteady effects of planet 
mass gain are of possible importance in that case. Consider 
the case of an initially steady accretion disk in which a planet 
of mass ratio q < lx 1CP 3 grows and increases its mass on a 
timescale shorter than the viscous timescale. We expect that 
the accretion rate onto the central star corresponds to the 
steady state rate at an earlier time when the planet mass was 
lower. Therefore, we expect that the accretion rate onto the 
star to be at least ~10% of the accretion flow rate outside the 
planet’s orbit in this mass range. 

For an inwardly migrating planet, the available time to 
produce an inner hole is limited by the planet’s migration 
timescale. Furthermore, a planet that moves with the accre¬ 
tion flow (Type II migration) would not be expected to have 
much effect on the accretion rate onto the central star. Cur¬ 
rent migration models in the planet mass range considered 
here (0.05 to 1 Mj), suggest that migration occurs on about 
the local viscous timescale (e.g., Bate et al. 2003). We then 



8 


Gas Flow Across Gaps in Protoplanetary Disks 



Fig. 2.— Plots of azimuthally-averaged density in units of E p versus radius in units of r p for the unevolved (dashed) and 
evolved (solid) models b (top-left), g (top-right), f (bottom-left), and 1 (bottom-right) in Table 1. 


expect that an inwardly migrating planet at some orbit ra¬ 
dius would have a more substantial interior disk and higher 
accretion rate onto the star than a planet on a fixed orbit at 
the same radius. Both the effects of planet mass gain and 
migration appear to lead to less depleted inner disks than the 
steady state models suggest. 


Present-day multi-dimensional hydrodynamical simula¬ 
tions are incapable of following a typical protostcllar disk 
to a steady state with spatial coverage down to the cen¬ 
tral star. To circumvent this problem, we have developed 
a one-parameter (namely the accretion efficiency) family of 
semi-analytic steady state solutions for the main disk, ex- 
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Density 



Density 



Fig. 3.— Density (in units described below) for a steady state disk containing an accreting planet, whose mass ratio is 
q = lx 10 —3 , is plotted as a function of radius in units of r p . The left plot (model b of Table 1) is for a constant disk viscosity 
and the right plot is for a viscosity that increases in radius as y/r (model f of Table 1). The bold solid line corresponds to a disk 
whose inner boundary lies on a star whose radius is very small compared to the orbital radius of the planet (r* = 0). The faint 
solid line is for an inner boundary r* = 0.4r p , typically applied in multidimensional hydrodynamical simulations. The dashed 
line is for the same disk without a planet. The unit of density in each plot is the density of the corresponding disk without a 
planet at r — r p . 


eluding the gap region. These solutions are combined with 
numerical simulations that include the gap region to provide 
approximate steady state global models. 

There is evidence in favor of the models being a reasonable 
representation of a steady state configuration. In particular, 
there is little indication of significant changes to the density 
profiles over the course of the numerical simulations (see Fig¬ 
ure 2). The density values at the inner boundary in the nu¬ 
merical simulations (which are free to change) change little. 
The steady state semi-analytic model makes a parameter-free 
prediction for the ratio of accretion onto the planet to the rate 
of accretion past the planet. The numerical simulations are 
in reasonable agreement with the predicted steady state ra¬ 
tio, generally within 25%. The effects of changing the radius 
of the star on the accretion rate obtained by the numerical 
models are also in accord with the predictions of the steady 
state analytic model. 

For a disk with thickness ratio H/r = 0.05 and planet to 
star mass ratios q of order 1 x 10~ 3 or less, about 10% or more 
of the accretion flow continues past the orbit of the planet 
and onto the central star. Such reductions would probably 
still typically produce detectable accretion onto the central 
star. Most of the accretion matter flows onto the planet. As 
mentioned earlier, we have generally assumed that the planet 
can accrete the material that reaches it, in accord with the 
expectations of current planet formation theory by Pollack et 
al. (1996) for the mass range we consider. If that assumption 
were incorrect to the extent that the planet does not accrete 
any mass, then the star should accrete at the rate expected 
in the absence of a q < 1 x 10 -3 planet. 

The effect of an accreting planet is to decrease the disk 
density over several times the planet’s orbital radius (see Fig¬ 
ure 3). The presence of an accreting planet can change the 


large-scale disk density distribution from a simple power law 
in radius. Inner circumstellar disks (interior to the orbit of 
the planet) occur in a steady state, even in the presence of a 
Jupiter-mass planet. The disk density inside the orbit of the 
planet (but outside the gap) is decreased as much as about 
an order of magnitude, due to the presence of a planet. Such 
inner disks may be missed or underestimated in numerical 
simulations whose inner boundaries are not sufficiently close 
to the star (see Figure 3). 

Producing cleaner inner holes with lower accretion rates 
onto the central star may be possible with mass ratios in 
excess of q = 1 x 10 -3 . For planet to star mass ratios q ~ 

5 x 10~ 3 , clean inner holes could result (Lubow, Seibert, & 
Artymowicz 1999). But, accretion may still occur at higher 
planet masses, if the disk becomes eccentric as a consequence 
of a tidal instability (Papaloizou, Nelson, & Masset 2001; Kley 

6 Dirksen 2005; D’Angelo, Lubow, & Bate 2005). 

For fixed disk viscosity and planet mass, a colder disk has 
less ability to penetrate the gap around a planet. The results 
(model h of Table 1) show that some reduction in flow through 
the gap occurs with a colder disk. For a mass ratio q = 
1 x 10~ 3 , the flow is reduced by about 25% when H/r is 
reduced from 0.05 to 0.03. 

In the case of the TW Hya, on the basis of low near-IR ex¬ 
cess, there is evidence of dust depletion inside 4 AU, although 
the depletion is not complete (Calvet et al. 2002). Based on 
emission, there is evidence of gas accretion that is one or 
two orders of magnitude lower than typical for much younger 
T Tauri stars (Muzerolle et al. 2000) . Some of the reduc¬ 
tion could be due to overall gas depletion associated with its 
~ 10 Myr age due to accretion. It is also possible that the 
presence of a planet of one Jupiter mass or less mass could 
play a role in decreasing, but not terminating, the accretion 
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in the inner region, by an order of magnitude. Similar con¬ 
siderations apply to the case of CoKu Tau/4, where there is 
evidence for dust depletion inside of 10 AU and no evidence of 
accretion (D’Alessio et al. 2005). We do not find that inner 
disk accretion can be suppressed by several orders of mag¬ 
nitude for the planet mass range considered (< lMj). Such 
suppression may be possible with a planet of several Mj. But 
as mentioned above, even this possibility is uncertain and re¬ 
quires further investigation. 
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